clear
I=imread(['case12_posta.jpg']);
imshow(I);
pause
prompt = {'Enter angle of rotation'};
dlgtitle = 'Input';
definput = {'0'};
dims = [1 40];
opts.Interpreter = 'tex';
answer = inputdlg(prompt,dlgtitle,dims,definput,opts);
angle=str2num(answer{1});
%angle=input('Input angle of rotation')
I = imrotate(I,angle);

BW=rgb2gray(I);
%threshold=[0.0375,0.0938];
%BW1=edge(BW,'canny',threshold,3);
%BW1=edge(BW,'canny')
[BW1,threshOut]=edge(BW,'canny');
imshow(BW1);
pause
BW1=edge(BW,'canny',threshOut,3);
imshow(BW1);
%pause
%imshow(BW2);
%yes=input('Filter looks ok? 0 for no, 1 for yes.');

%Do a whole if loop to check that they are ok.
pause
BW2=bwareaopen(BW1,20);
imshow(BW1);

%BW2=BW1;
imshow(BW2);
pause
f = msgbox("Program is paused. Zoom into thh bifurcation and press any key.");
%prompt=('Program is paused. Zoom into the bifurcation to make ROI easier then press spacebar.')
pause
%prompt=('Draw a line above the main vessel to remove unwanted objects')
%[xtop,ytop]=getpts
xlabel('Eliminate unwanted noise around edges with ROI');
h = drawpolygon(gca,'FaceAlpha',0);
h.Color = 'yellow';
h2=drawpolygon(gca,'FaceAlpha',0);
h2.Color= 'yellow';
h3=drawpolygon(gca,'FaceAlpha',0);
h3.Color='yellow';

mask=createMask(h,BW2);
mask2=createMask(h2,BW2);
mask3=createMask(h3,BW2);

BW2(mask)=[0];
BW2(mask2)=[0];
BW2(mask3)=[0];

%BW2=bwareaopen(BW2,20);
imshow(BW2)
f = msgbox("Program is paused, press spacebar to continue");
%prompt=('Program is paused, press spacebar to continue')
pause
xlabel('Draw box around bifurcation area, placing 4 points with the mouse')
%f = msgbox("Draw a box around the bifurcation area, placing 4 points with the mouse");
%prompt=('Draw a box around the bifurcation area, placing 4 points with the mouse');
%Draw box around bifurcation, starting lower left hand corner and going
%clockwise

[xk yk]=getpts
imshow(I);
hold on

%Get the box ready
maxx=(max(xk));
maxx=round(maxx);
minx=round(min(xk));
maxy=round(max(yk));
miny=round(min(yk));
ysize=(maxx-minx);
xsize=(maxy-miny);

[numRows,numColumns]=size(BW2)

 for i=miny:1:maxy 
  for j=minx:1:maxx
      if BW2(i,j)==1
          plot(j,i,'*','MarkerSize',1,'Color',[99/255,141/255,187/255]);
      end
      j=j+1;
  end
  i=i+1;
 end

xlabel('Program paused. Press any key to continue')
pause 

%Select Points around ostium
%prompt=('Select 3 points along vessel edge, 1 on MV, 1 approximately on ostium, 1 in SB.')
%[linex,liney]=getpts


%[line1x,line1y]=pointextract2(linex(1),liney(1),linex(2),liney(2),BW2);
%[skelatonx2,skelatony2]=pointextract2(linex(2),liney(2),linex(3),liney(3),BW2);
%{
x=[linex(2),linex(3)];
y=[liney(2),liney(3)];
poly1=polyfit(x,y,1);
x1=round(x(1));
x2=round(x(2));
y1=round(y(1));
y2=round(y(2));
skelatonx2=[1];
skelatony2=[1];
for j=y2:1:y1
        checkpoint=round((j-poly1(2))/poly1(1))-4;
        for i=1:1:8
            if BW2(j,checkpoint)==1
                skelatonx2=[skelatonx2,checkpoint];
                skelatony2=[skelatony2,j];
            end
            i=i+1;
            checkpoint=checkpoint+1;
        end
        j=j+1;
end
%}
%delete duplicates
%yo=size(line1x);
%m=size(skelatonx2);
%if skelatonx2(10)<skelatonx2(1)
%skelatonx2=flip(skelatonx2);
%skelatony2=flip(skelatony2);
%end
%if line1x(10)<line1x(1)
%    line1x=flip(line1x);
%    line1y=flip(line1y);
%end
%for i=1:1:8
%    checkx=skelatonx2(i);
%    checky=skelatony2(i);
%    for j=1:1:yo(1,2)
%        max1=line1x(j);
%        may1=line1y(j);
%        if checkx==max1 && checky==may1
%            skelatonx2=skelatonx2((i+1):end);
%            skelatony2=skelatony2((i+1):end);%
%end
%        j=j+1;
%    end
%    i=i+1;
%end
%plot(line1x,line1y,'o','Color',[56/255,124/255,68/255]);
%plot(skelatonx2,skelatony2,'*','Color',[56/255,124/255,68/255]);

%concatenatex=[line1x,skelatonx2];
%concatenatey=[line1y,skelatony2];

%slope=[1];
%heifer=size(concatenatex);
%heifer2=heifer(1,2);
%heifer3=heifer2-12;
%for i=1:1:heifer3
%    value=i+10;
%    shalomx=concatenatex(i:1:value);
%    shalomy=concatenatey(i:1:value);
%    sloop=polyfit(shalomx,shalomy,1);
%    slope=[slope,sloop(1)];
%    i=i+1;
%end
%slope=slope(2:end)

%slopeMV=(slope(1)+slope(2)+slope(3)+slope(4)+slope(5))/5;
%slopeSB=(slope(heifer3)+slope(heifer3-1)+slope(heifer3-2)+slope(heifer3-3)+slope(heifer3-4))/5;
%avgslope=(slopeMV+slopeSB)/2;
%if slope(1)>avgslope
%    for i=1:1:heifer3
%        comp=slope(i);
%        if abs(comp)>abs(avgslope)
%            amigo=i;
%            break
%        end
%    end
%else
%for i=1:1:heifer3
%    comp=slope(i);
%    if comp>avgslope
%        amigo=i;
%        break
%    end
%    i=i+1;
%end
%end
%%%
%amigo=amigo+5;
%plot(concatenatex(1),concatenatey(1),'o','Color',[56/255,124/255,68/255]);
%plot(concatenatex(amigo),concatenatey(amigo),'o','Color',[152/255,1,152/255]);
%checkostium=input('Is the ostium location ok? 0 for no, 1 for yes.')
%if checkostium==0
%    amigo=amigo-5;
%    plot(concatenatex(amigo),concatenatey(amigo),'x','Color',[152/255,1,152/255]);
%    checkostium=input('Is the new ostium ok? (x)');
%    if checkostium==0
%        amigo=amigo+10;
%        plot(concatenatex(amigo),concatenatey(amigo),'x','Color',[152/255,1,152/255]);
%        checkostium=input('Is the new ostium (x) ok?');
%        if checkostium==0
xlabel('Select side branch ostium point')
            prompt=('Select Ostium point manually.');
            %[concatenatex(amigo),concatenatey(amigo)]=getpts;
            [ostiumx,ostiumy]=getpts;
            %plot(concatenatex(amigo),concatenatey(amigo),'o','Color',[152/255,1,152/255]);
            plot(ostiumx,ostiumy,'o','Color',[152/255,1,152/255]);
%        end
%   end
%end

%%
%Centerline of MV
xlabel('Select two points for upper and lower main vessel, along a straight section, undiseased section')
[MVx MVy]=getpts;
MVtopx1=MVx(1);
MVtopx2=MVx(2);
MVtopy1=MVy(1);
MVtopy2=MVy(2);

MVbotx1=MVx(3);
MVbotx2=MVx(4);
MVboty1=MVy(3);
MVboty2=MVy(4);

[skelMVtopx, skelMVtopy]=pointextract2(MVtopx1,MVtopy1,MVtopx2,MVtopy2,BW2);
[skelMVbotx, skelMVboty]=pointextract2(MVbotx1,MVboty1,MVbotx2,MVboty2,BW2);
skelMVtop=polyfit(skelMVtopx,skelMVtopy,1);
skelMVbot=polyfit(skelMVbotx,skelMVboty,1);
plotMVtopx=linspace(MVtopx1,MVtopx2);
plotMVbotx=linspace(MVbotx1,MVbotx2);
plot(plotMVtopx,polyval(skelMVtop,plotMVtopx),'Color',[191/255,209/255,229/255]);
plot(plotMVbotx,polyval(skelMVbot,plotMVbotx),'Color',[191/255,209/255,229/255]);
plot(skelMVtopx,skelMVtopy,'*','Color',[6/255,77/255,135/255]);
plot(skelMVbotx,skelMVboty,'*','Color',[6/255,77/255,135/255]);

%Estimate Centerline
[centxy]=lineintersect(skelMVtop(1),skelMVtop(2),skelMVbot(1),skelMVbot(2));
mcent=(skelMVtop(1)+skelMVbot(1))/2;
[centeqn]=pointslope(centxy(1),centxy(2),mcent);
plot(plotMVtopx,polyval(centeqn,plotMVtopx),'Color',[191/255,209/255,229/255]);
%%
%Centerline of MB
xlabel('Select two points for upper & lower main main branch; outside points first. (nondiseased section)')
[MBx MBy]=getpts;
MBtopx1=MBx(1);
MBtopx2=MBx(2);
MBtopy1=MBy(1);
MBtopy2=MBy(2);

MBbotx1=MBx(3);
MBbotx2=MBx(4);
MBboty1=MBy(3);
MBboty2=MBy(4);

[skelMBtopx, skelMBtopy]=pointextract2(MBtopx1,MBtopy1,MBtopx2,MBtopy2,BW2);
[skelMBbotx, skelMBboty]=pointextract2(MBbotx1,MBboty1,MBbotx2,MBboty2,BW2);
skelMBtop=polyfit(skelMBtopx,skelMBtopy,1);
skelMBbot=polyfit(skelMBbotx,skelMBboty,1);
plotMBtopx=linspace(MBtopx1,MBtopx2);
plotMBbotx=linspace(MBbotx1,MBbotx2);
plot(plotMBtopx,polyval(skelMBtop,plotMBtopx),'Color',[191/255,209/255,229/255]);
plot(plotMBbotx,polyval(skelMBbot,plotMBbotx),'Color',[191/255,209/255,229/255]);
plot(skelMBtopx,skelMBtopy,'*','Color',[6/255,77/255,135/255]);
plot(skelMBbotx,skelMBboty,'*','Color',[6/255,77/255,135/255]);

%Estimate Centerline
[centMB]=lineintersect(skelMBtop(1),skelMBtop(2),skelMBbot(1),skelMBbot(2));
mcentMB=(skelMBtop(1)+skelMBbot(1))/2;
[centeqnMB]=pointslope(centMB(1),centMB(2),mcentMB);
plot(plotMBtopx,polyval(centeqnMB,plotMBtopx),'Color',[191/255,209/255,229/255]);
%%
%Centerline of SB
xlabel('Select two points for upper and lower side branch; select outside points first, nondiseased portion')
[SBx SBy]=getpts;
SBtopx1=SBx(1);
SBtopx2=SBx(2);
SBtopy1=SBy(1);
SBtopy2=SBy(2);

SBbotx1=SBx(3);
SBbotx2=SBx(4);
SBboty1=SBy(3);
SBboty2=SBy(4);

[skelSBtopx, skelSBtopy]=pointextract2(SBtopx1,SBtopy1,SBtopx2,SBtopy2,BW2);
[skelSBbotx, skelSBboty]=pointextract2(SBbotx1,SBboty1,SBbotx2,SBboty2,BW2);
skelSBtop=polyfit(skelSBtopx,skelSBtopy,1);
skelSBbot=polyfit(skelSBbotx,skelSBboty,1);
plotSBtopx=linspace(SBtopx1,SBtopx2);
plotSBbotx=linspace(SBbotx1,SBbotx2);
plot(plotSBtopx,polyval(skelSBtop,plotSBtopx),'Color',[191/255,209/255,229/255]);
plot(plotSBbotx,polyval(skelSBbot,plotSBbotx),'Color',[191/255,209/255,229/255]);
plot(skelSBtopx,skelSBtopy,'*','Color',[6/255,77/255,135/255]);
plot(skelSBbotx,skelSBboty,'*','Color',[6/255,77/255,135/255]);

%Estimate Centerline
[centSB]=lineintersect(skelSBtop(1),skelSBtop(2),skelSBbot(1),skelSBbot(2));
mcentSB=(skelSBtop(1)+skelSBbot(1))/2;
[centeqnSB]=pointslope(centSB(1),centSB(2),mcentSB);
plot(plotSBtopx,polyval(centeqnSB,plotSBtopx),'Color',[191/255,209/255,229/255]);

%% Diameter Calculation
[MVpoints] = finddiameter(skelMVtop,skelMVbot,centeqn,skelMVtopx,skelMVtopy,skelMVbotx,skelMVboty);
[MBpoints] = finddiameter(skelMBtop,skelMBbot,centeqnMB,skelMBtopx,skelMBtopy,skelMBbotx,skelMBboty);
[SBpoints] = finddiameter(skelSBtop,skelSBbot,centeqnSB,skelSBtopx,skelSBtopy,skelSBbotx,skelSBboty);

[MVDiameter] = finddistance(MVpoints)
[MBDiameter] = finddistance(MBpoints)
[SBDiameter] = finddistance(SBpoints)
%% Angles !!
xlabel('Draw MV centerline')
billybob = drawline();
pause
MVL = billybob.Position
delete(billybob)
xMVpt=[MVL(1,1),MVL(2,1)];
yMVpt=[MVL(1,2),MVL(2,2)];
xMV = linspace(MVL(1,1),MVL(2,1));
MVform = polyfit(xMVpt,yMVpt,1)
plot(xMV,polyval(MVform,xMV),'--','Color','b');

xlabel('Draw MB centerline')
billybob = drawline();
pause
MBL = billybob.Position
delete(billybob)
xMBpt=[MBL(1,1),MBL(2,1)];
yMBpt=[MBL(1,2),MBL(2,2)];
xMB = linspace(MBL(1,1),MBL(2,1));
MBform = polyfit(xMBpt,yMBpt,1)
plot(xMB,polyval(MBform,xMB),'--','Color','b')

xlabel('Draw SB centerline')
billybob = drawline();
pause
SBL = billybob.Position
delete(billybob)
xSBpt=[SBL(1,1),SBL(2,1)];
ySBpt=[SBL(1,2),SBL(2,2)];
xSB = linspace(SBL(1,1),SBL(2,1));
SBform = polyfit(xSBpt,ySBpt,1)
plot(xSB,polyval(SBform,xSB),'--','Color','b')

%% Carina Point
xlabel('Select carina point')
[carx,cary]=getpts;
%Define carina Line
[carinaline]=pointslope(carx,cary,-1/MVform(1));
%First guess for carina diameter
[car1]=lineintersect(carinaline(1),carinaline(2),SBform(1),SBform(2));
[car2]=lineintersect(carinaline(1),carinaline(2),MVform(1),MVform(2));
%Loops for getting exact carina diameter intersection
car1xa=round(car1(1));
car1ya=round(car1(2));
car1xb=round(car1(1));
car1yb=round(car1(2));

car2xa=round(car2(1));
car2ya=round(car2(2));
car2xb=round(car2(1));
car2yb=round(car2(2));
j=1
for i=1:1:10
    if car1ya<0
        car1ya=-car1ya;
    end
    if car1xa<0
        car1xa=-car1xa;
    end
   check=BW2(car1ya,car1xa)
   if check==1
       car1=[car1xa,car1ya];
       break
   end
   check2=BW2(car1yb,car1xb)
   if check2==1
       car1=[car1xb,car1yb]
       break
   end
   tester=(j)/sqrt(1+(carinaline(1))^2);
   test1x=car1(1)+tester;
   test2x=car1(1)-tester;
   test1y=polyval(carinaline,test1x);
   test2y=polyval(carinaline,test2x);
   car1xa=round(test1x);
   car1ya=round(test1y);
   car1xb=round(test2x);
   car1yb=round(test2y);
   j=j+1;
   if test1x>514
       test1x=0;
   end
   if test1y>514
       test1y=0;
   end
   if test2x>514
       test2x=0;
   end
   if test2y>514
       test2y=0;
   end
   plot(test1x,test1y,'x','Color',[0.5 0 0.8]);
   plot(test2x,test2y,'x','Color',[0.5 0 0.8]);
  %plot(car1xa,car1ya,'o','Color',[0.5 0 0.8]);
   i=i+1;
end
for i=1:1:10
   check=BW2(car2ya,car2xa)
   if check==1
       car2=[car2xa,car2ya];
       break
   end
   check2=BW2(car2yb,car2xb)
   if check2==1
       car2=[car2xb,car2yb];
       break
   end
   tester=(i)/sqrt(1+(carinaline(1))^2);
   test1x=car2(1)+tester;
   test2x=car2(1)-tester;
   test1y=polyval(carinaline,test1x);
   test2y=polyval(carinaline,test2x);
   car2xa=round(test1x);
   car2ya=round(test1y);
   car2xb=round(test2x);
   car2yb=round(test2y);
   i=i+1;
   plot(test1x,test1y,'o','Color',[0.5 0 0.8]);
   plot(test2x,test2y,'o','Color',[0.5 0 0.8]);
end
%plot actual carina lines
carinalinex=linspace(car2(1),car1(1));
plot(carinalinex,polyval(carinaline,carinalinex),'Color',[130/255,0,86/255]);
plot(car2(1),car2(2),'d','Color',[228/255,0,150/255]);
plot(car1(1),car1(2),'d','Color',[228/255,0,150/255]);

%allow user intervention
 choice = questdlg('Carina diameter ok?','Choose','Yes','No','No');
    if isequal(choice,'No')
        xlabel('Select two points where carina line intersects vessel edges')
        [cardx,cardy]=getpts;
        car2(1)=cardx(1);
        car2(2)=cardy(1);
        car1(1)=cardx(2);
        car1(2)=cardy(2);
        plot(car2(1),car2(2),'d','Color',[1,43/255,0])
        plot(car1(1),car1(2),'d','Color',[1,43/255,0])
    end
%{
confirm=input('Carina diameter ok? 1=yes, 0=no');
if confirm==0
    xlabel('Select Carina points manually');
    [cardx,cardy]=getpts;
    car2(1)=cardx(1);
    car2(2)=cardy(1);
    car1(1)=cardx(2);
    car1(2)=cardy(2);
    plot(car2(1),car2(2),'d','Color',[1,43/255,0])
    plot(car1(1),car1(2),'d','Color',[1,43/255,0])
end
%}
CarinaLength=pointdistance(car2(1),car2(2),car1(1),car1(2))
%Find points where carina diameter line intersects with the outside lines
%Display points, allow manual adjustment
%% Ostium Distance
%Ok, use ostium point and perpendicular slope to get line for ostium
%location, and highlight ostium line intersection with MV centerline. Then
%find distance between ostium line and carina diameter line.
[ostiumline]=pointslope(ostiumx,ostiumy,carinaline(1));
%[ostiumline]=pointslope(concatenatex(amigo),concatenatey(amigo),carinaline(1));
[xdistancepoint]=lineintersect(ostiumline(1),ostiumline(2),centeqn(1),centeqn(2));
%xostiumline=linspace(concatenatex(amigo),xdistancepoint(1));
xostiumline=linspace(ostiumx,xdistancepoint(1));
plot(xostiumline,polyval(ostiumline,xostiumline),'--','Color',[56/255,124/255,68/255]);
plot(xdistancepoint(1),xdistancepoint(2),'d','Color',[152/255,1,152/255]);
plot(carx,cary,'d','Color',[228/255,0,150/255]);
XLength=linedistance(carinaline(1),carinaline(2),ostiumline(1),ostiumline(2))
%% Diameter Calculation
[MVpoints] = finddiameter(skelMVtop,skelMVbot,centeqn,skelMVtopx,skelMVtopy,skelMVbotx,skelMVboty);
[MBpoints] = finddiameter(skelMBtop,skelMBbot,centeqnMB,skelMBtopx,skelMBtopy,skelMBbotx,skelMBboty);
[SBpoints] = finddiameter(skelSBtop,skelSBbot,centeqnSB,skelSBtopx,skelSBtopy,skelSBbotx,skelSBboty);

[MVDiameter] = finddistance(MVpoints)
[MBDiameter] = finddistance(MBpoints)
[SBDiameter] = finddistance(SBpoints)


%% Angle calculation based on centerlines
%Calculate Angles Between Branches:
%A=angle between main branch and side branch
%B=angle between main vessel and main branch
%C=angle between main vessel and side branch
A=atand((MBform(1)-SBform(1))/(1+MBform(1)*SBform(1)))
B=atand((MBform(1)-MVform(1))/(1+MBform(1)*MVform(1)))
C=atand((MVform(1)-SBform(1))/(1+MVform(1)*SBform(1)))

xlabel('Program paused. Press any key to continue')
pause
%% Catheter Diameter
hold off
imshow(I);
xlabel('Select 4 points around the catheter region')
prompt=('Select 4 points around the catheter')
[xk yk]=getpts
imshow(I);
hold on

%Get the box ready
maxx=(max(xk));
maxx=round(maxx);
minx=round(min(xk));
maxy=round(max(yk));
miny=round(min(yk));
ysize=(maxx-minx);
xsize=(maxy-miny);

[numRows,numColumns]=size(BW2)

 for i=miny:1:maxy 
  for j=minx:1:maxx
      if BW2(i,j)==1
          plot(j,i,'b*','MarkerSize',1);
      end
      j=j+1;
  end
  i=i+1;
 end
xlabel('Program paused. Press any key to continue');
pause 
xlabel('Select 2 points on edge of the catheter, and two points on another edge.');
prompt=('Select 2 points on one edge of the catheter, and two points on another edge.');
[cathx,cathy]=getpts


[cath1x,cath1y]=pointextract2(cathx(1),cathy(1),cathx(2),cathy(2),BW2);
[cath2x,cath2y]=pointextract2(cathx(3),cathy(3),cathx(4),cathy(4),BW2);
cath1=polyfit(cath1x,cath1y,1);
cath2=polyfit(cath2x,cath2y,1);

[cathcentpoint]=lineintersect(cath1(1),cath1(2),cath2(1),cath2(2));
mcath=(cath1(1)+cath2(1))/2;
[cathcenteqn]=pointslope(cathcentpoint(1),cathcentpoint(2),mcath);

[cathpoints] = finddiameter(cath1,cath2,cathcenteqn,cath1x,cath1y,cath2x,cath2y);
[cathpixels] = finddistance(cathpoints);

%Draw Lines around diameter
halfdiameter=cathpixels/2;
cathperpslope=-1/cathcenteqn(1);
halfcath=(halfdiameter)/sqrt(1+(cathperpslope)^2);
dummyx=(cathx(1)+cathx(4))/2;
dummyy=polyval(cathcenteqn,dummyx);
cathperp=pointslope(dummyx,dummyy,cathperpslope);

xplus=dummyx+halfcath;
yplus=polyval(cathperp,xplus);

xminus=dummyx-halfcath;
yminus=polyval(cathperp,xminus);

pluseqn=pointslope(xplus,yplus,cathcenteqn(1));
minuseqn=pointslope(xminus,yminus,cathcenteqn(1));
xplotplus=linspace(xplus-8,xplus+8);
xplotminus=linspace(xminus-8,xplus+8);

plot(xplotplus,polyval(pluseqn,xplotplus),'y');
plot(xplotminus,polyval(minuseqn,xplotminus),'y');

%Get French Size and convert diameter values
prompt = {'Enter Catheter French size (number only)'};
dlgtitle = 'Input';
definput = {'5'};
dims = [1 40];
opts.Interpreter = 'tex';
answer = inputdlg(prompt,dlgtitle,dims,definput,opts);
frenchsize=str2num(answer{1});
%frenchsize=input('Input Catheter French Size (number only)');
multiplier=frenchsize/cathpixels;
FrMVDiameter=MVDiameter*multiplier
FrMBDiameter=MBDiameter*multiplier
FrSBDiameter=SBDiameter*multiplier
FrCarinaLength=CarinaLength*multiplier
FrXLength=XLength*multiplier

%Convert from French Size to mm
mmperfrench=1/3;
MVDiameterMM=FrMVDiameter*mmperfrench;
MBDiameterMM=FrMBDiameter*mmperfrench;
SBDiameterMM=FrSBDiameter*mmperfrench;
CarinaLengthMM=FrCarinaLength*mmperfrench;
XLengthMM=FrXLength*mmperfrench;
%% Output Diameters,angles, arrays, and point matrices to Excel spreadsheet
%so that results can be 'recorded' and replotted easily if needed.

labels=["Angle A";"Angle B"; "Angle C"; "MV DIAM"; "MBDIAM";"SBDIAM";
    "CarinaLength";"X-Length";"CathWidth";"MVDIAM (adj";"MBDIAM (adj)";
    "SBDIAM (adj)";"CARDIAM (adj)";"XLength (Adj)";"Cath Width adj";"Cath Size (french"
    "MVDIAM (MM)"; "MBDIAM (MM)";"SBDIAM (MM)";"CARDIAM (MM)";"XLength (MM)"];
values=[A;B;C;MVDiameter;MBDiameter;SBDiameter;CarinaLength;XLength;
    cathpixels;FrMVDiameter;FrMBDiameter;FrSBDiameter;FrCarinaLength;
    FrXLength;frenchsize;frenchsize;MVDiameterMM;MBDiameterMM;SBDiameterMM;CarinaLengthMM;XLengthMM];
table(labels,values)
pause
table(values)

spy
xlabel('Cheers')
%choice = questdlg('Did you have fun','Choose','Yes','No','No');
 %   if isequal(choice,'Yes')
  %      spy
   %     xlabel('Good')
    %else
     %   why
    %end
%%
function [skelatonx, skelatony]=pointextract2(x1,y1,x2,y2,pointmap)
x=[x1 x2];
y=[y1 y2];
poly1=polyfit(x,y,1);
x1=round(x1);
x2=round(x2);
y1=round(y1);
y2=round(y2);
if x1>x2
    placehold=x1;
    x1=x2;
    x2=placehold;
    placehold=y1;
    y1=y2;
    y2=placehold;
end

skelatonx=[1];
skelatony=[1];

m=abs(poly1(1));
checkx=x1;
if m<1
for j=x1:1:x2
    checkpoint=round(polyval(poly1,checkx)-5);
    for i=1:1:10
        if pointmap(checkpoint,checkx)==1
            skelatonx=[skelatonx checkx];
            skelatony=[skelatony checkpoint];
        end
        i=i+1;
        checkpoint=checkpoint+1;
    end
    j=j+1;
    checkx=checkx+1;
end
else
    if y1>y2
    placehold=x1;
    x1=x2;
    x2=placehold;
    placehold=y1;
    y1=y2;
    y2=placehold;
    end
  for j=y1:1:y2
      checkpoint=round((j-poly1(2))/poly1(1))-4;
      for i=1:1:8
          if pointmap(j,checkpoint)==1
              skelatonx=[skelatonx,checkpoint];
              skelatony=[skelatony,j];
          end
          i=i+j;
          checkpoint=checkpoint+1;
      end
      j=j+1;
  end
end
%pointfit=polyfit(skelaton(2:end),skelatony(2:end),1);
skelatonx=skelatonx(2:end);
skelatony=skelatony(2:end);
end

function [lineeqn]=pointslope(x1,y1,m)
b=y1-m*x1;
lineeqn=[m,b]; 
end


%Function: line intersect
%Input: m1,b1,m2,b2
%Output: x,y point where the two lines intersect
function [xy]=lineintersect(m1,b1,m2,b2)
x=(b2-b1)/(m1-m2);
y=m1*x+b1;
xy=[x,y];
end
function[distance]=pointdistance(x1,y1,x2,y2)
distance=sqrt((x2-x1)^2+(y2-y1)^2);
end

%% How To Use finddiameter!
%Inputs: (top equation, bottom equation, middle equation, top array (x and
%y), bottom array (x and y))
%Outputs: An array with the following format:
%     x1top y1top x1bot y1bot
%     x2top y2top x2bot y2bot
%     x3top y3top x3bot y3bot
%Then you should average the three distances between the three pairs of points. ! :)
function [points] = finddiameter(toplineequation,bottomlineequation,skelatonlineequation,toparrayx,toparrayy,bottomarrayx,bottomarrayy)
bskel=skelatonlineequation(2);
mskel=skelatonlineequation(1);
mperp=-1/mskel;
btop=toplineequation(2);
mtop=toplineequation(1);
bbot=bottomlineequation(2);
mbot=bottomlineequation(1);
%x1
top1x=toparrayx(2);
top1y=toparrayy(2);
bot1x=bottomarrayx(2);
bot1y=bottomarrayy(2);
btop1=top1y-mperp*top1x;
bbot1=bot1y-mperp*bot1x;
xtop1=(btop1-bskel)/(mskel-mperp);
xbot1=(bbot1-bskel)/(mskel-mperp);
if xtop1>xbot1
    x1=xtop1;
else
    x1=xbot1;
end
%x2
top2x=toparrayx(end);
top2y=toparrayy(end);
bot2x=bottomarrayx(end);
bot2y=bottomarrayy(end);
btop2=top2y-(mperp*top2x);
bbot2=bot2y-(mperp*bot2x);
xtop2=(btop2-bskel)/(mskel-mperp);
xbot2=(bbot2-bskel)/(mskel-mperp);
if xtop2<xbot2
    x2=xtop2;
else
    x2=xbot2;
end
%x3
x3=(x1+x2)/2;

%y1
y1=mskel*x1+bskel;
bp1=y1-mperp*x1;
x1top=(btop-bp1)/(mperp-mtop);
y1top=mtop*x1top+btop;
x1bot=(bbot-bp1)/(mperp-mbot);
y1bot=mbot*x1bot+bbot;

%y2
y2=mskel*x2+bskel;
bp2=y2-mperp*x2;
x2top=(btop-bp2)/(mperp-mtop);
y2top=mtop*x2top+btop;
x2bot=(bbot-bp2)/(mperp-mbot);
y2bot=mbot*x2bot+bbot;

%y3
y3=mskel*x3+bskel;
bp3=y3-mperp*x3;
x3top=(btop-bp3)/(mperp-mtop);
y3top=mtop*x3top+btop;
x3bot=(bbot-bp3)/(mperp-mbot);
y3bot=mbot*x3bot+bbot;

points=[x1top y1top x1bot y1bot
    x2top y2top x2bot y2bot
    x3top y3top x3bot y3bot];
end

function [distance] = finddistance(diameterarray)
x1a=diameterarray(1,1);
y1a=diameterarray(1,2);
x1b=diameterarray(1,3);
y1b=diameterarray(1,4);

x2a=diameterarray(2,1);
y2a=diameterarray(2,2);
x2b=diameterarray(2,3);
y2b=diameterarray(2,4);

x3a=diameterarray(3,1);
y3a=diameterarray(3,2);
x3b=diameterarray(3,3);
y3b=diameterarray(3,4);

dist1=sqrt(((x1a-x1b)^2)+((y1a-y1b)^2));
dist2=sqrt(((x2a-x2b)^2)+((y2a-y2b)^2));
dist3=sqrt(((x3a-x3b)^2)+((y3a-y3b)^2));

distance=(dist1+dist2+dist3)/3;
end

function [distance]=linedistance(m1,b1,m2,b2)
%Find the distance between two parallel lines
distance=(abs(b1-b2))/(sqrt(1+(m1^2)));
end
